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SUMMARY 


For steady state time harmonic problems (rigid format 8), the NASTRAN pro- 
gram as currently configured treats internal structural damping through the 
introduction of a singlo structural element damping coefficient that typically 
is viewed as the ratio of the complex to real modulus of elasticity (E^/E^) , 

For problems dealing with two or three dimensional dynamic linear viscoelastic- 
ity (fe.g. a Kelvin-Voigc viscoelastic model), the present NASTRAN capability 
cannot directly handle this situation wherein two independent damping coeffici- 
ents are required to properly model the dissipation phenomenon. A technique is 
presented whereby the user can adapt the standard versions of NASTRAN (without 
resorting to either DMAP and/or FORTRAN coding changes) for the purpose of 
treating this class of problem. 


INTRODUCTION 


This paper is concerned with the solution to ] , 2, or 3 dimensional steady 
state (time harmonic) structural response problems wherein part or all of the 
structure is comprised of a linear viscoelastic material. In particular, 
attention is focused on the representation of the viscoelastic dissipation (or 
equivalently sound absorption) properties of this class of materials. Typi- 
cally, rubber-like materials fall into the category of interest. In situations 
where the driving frequencies of the applied loading is large, the effect of 
the energy dissipation characteristics on the overall dynamic response (partic- 
ularly in wave propagation problems) can be significant. Consequently, it is 
important that the material physical properties are modeled as accurately as 
possible. In this paper, the Kelvln-Voigt viscoelastic model is selected 
wherein the corresponding continuous field equations are given by (ref. (1)) 
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where X^,p are the usual Lame" elastic constants; X^, are their viscous 
counterparts; p denotes the mass per unit volume, u is the displacement vector; 
t is time and V is the vector "del" operator (3/Sx, 3/ay, 3/3z). The steady 
state harmonic version of equation (1) is obtained by substituting 
“ — ’i’XCxit 

u = u^e into equation (1) which results in the form 
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where and are the complex wave speeds given by the expressions 


c2 = c2(l+in^) 
c2 a c2(l+irij) 


(3) 


In equations (3) , c and c . are the usual elastic shear and dilatatlonal wave 

s cl 

speeds, and and are the associated sherr and dilauational dissipation 

constants; these four constants are related to the constants originally employed 
in equation (1) through the relations; 
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where X and . are two independent frequency dependent viscoelastic constants 
which are determined experimentally. The X^, constants are related to the 


“i 'i i 

X , y constants of Eq (1) through the relations X 


(1) 


X^ and yi = w y^. 
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Depending on one's viewpoint, the pair of parameter X , y' (or alternatively 
n » Hj) can be viewed as the two independent parameters which describe the 

S Q 

dissipation characteristics of the viscoelastic media. Further, the independ- 
ent constants A , (or alternatively c^, c^) can be viewed as the two 
independent constants which define the elastic characteristics of the media. 


The current version of NASTRAN can treat a solid media governed by 
equation (2) by using 2 and 3 dimensional elements (e.g. CTRMEM, CQDMEM, 

CTRAPGR, CTETRA, CTRIARG to name a few) in conjunction with both rigid format 8 
and the introduction of a loss factor (e.g. the GE input variable on a MATl or 
MAT2 card; this is referred to as the '’structural element imping coefficient” 
in the NASTRAN manual). Unfortunately, however, the user can introduce only one 
independent loss factor. A single, rather than two, independent loss factor can 
properly represent the Kelvin-Voigt model if the relation 



is met. Substituting equation (6) 
fact that the Young's modulus, E, 
yields the relations 


into equations (5) in conjunction with the 
is related to p, A through E ™ p (3X+2p) / (X+p) 

rig = eVe"" (7) 
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Thus, the more limited sir.gle parameter loss factor case can directly be 
implenientated in NASTRAN by assuming equation (6) holds and setting GE = on 
a MAT 1 or MAT 2 card as appropriate. 

The remainder of this paper is concerned with the situation wherein the 
viscoelastic physical -onstants are such that n, q (l.e. the special case 

G S 

implied by equations (6) are not satisfied). A technique is presented which 
permits the user to treat this more general case without having to resort to 
either DMAP instructions and/or modified FORTRAN coding modifications to the 
original NASTRAN program. Because of this general type "fix", the approach 
also has applications to other finite element programs having the same one 
parameter type dissipation limitation found in NASTRAN. 


IMPLEMENTATION OF TWO PARAMETER LOSS FACTORS 


Here we consider finite elements representing equations (2) for the case 
n n I • The finite element formulation representing equations (2) leads to a 

complex set of simultaneous equations to be solved of the form (ref. (2)) 
[K]{U}={P} (8) 

where [K]-[-u)^^M]+ito[B] + [K] ] (9) 


where {U} is the vector of nodal displacement amplitudes, {P} Is the vector of 
applied forces, [M] , is the assembled mass matrix, [B] is the assembled damping 
matrix and [K] is the elastic stiffness matrix. The formation of the mass, 
stiffness and non-structural damping portion of the matrices in Eq. (8) fol ow 
exactly the same process used in modeling elasticity type problems, hence will 
not be commented on here (e»g* see refs. (2,3). Attention is consequently 
focused on the formation of the structural damping part of the [K] matrix; let 


us define as the individual element structural damping portion due to the 

contiibution of the e^^'^ element- NASTRAN currently forms [K]^ from the 
relation ^ 

[C]"^njj[G]^[c;Jdxdydz (10) 


[K] 


e e 
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where [K] is the individual element elastic stiffness matrix, [C] is the 

^ r 

corresponding displacement-strain matrix and [G] is the elastic stress-strain 

r ^ 

law matrix. For isotropic materials, is of the form; 
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[G]' 


where a 


11 ’ 22 


, . . and e 


11 ’ 22 


... are the element stresses and strains. 


So far we have described, via equations (10) and (11), the form currently 
in NASTRAN. Next we consider the form of [K]^ we would like to have in order 

to implement the two independent parameter formulation. By comparing what we 
have to what we would like to have, the ’’fix" to NASTRAN will become evident. 

We start by making the important observation that the desired general continu- 
ous viscoelastic form of equation (2) can be derived from the usual elastic 
derivation for the dynamic equations of elasticity (ref. (4)). This is 
accomplished by following the elastic derivation of ref. (4) with the modifica- 
tion that, in place of the usual isotropic stress-strain law (i.e. Eq. (11)), 
we use 


{a}=[[G]’^ + i[G]^]{s} 
e e 


( 12 ) 


iL V 

where [G]^ is the same expression as [G]^ except all r superscripts are replaced 

with i superscripts for the A, p entries* Similarly, in the finite element 

formulation, all we need do is replace the usual [G]^ matrix in the stiffness 

derivation with [G] ^+i [G] 

e e 


Thus, the element stiffness becomes 


[K], = 
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[C]'^[[G]^+i[G]^][C] dxdyd z 
e e 
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[C]"^[G]^[C]dxdydz + i 
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[C ]^[G]^[C]dxdydz 


(12b) 


Term 1. . .Usual elastic stiffness 
contribution to [K] of Eq. (8). 


Term 2 .. .viscoelastic dissipation 
contribution to [K] of equation (8) 
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Nt»(f w 0 are in a poa'tion uo build the desired two .Independent parameter visco- 
elastic finite element. Tlie total complex element contribution, to the 

assembled [K] matrix In equation (9) is formed in two parts through the creation 
of an overlapping double viscoelastic element. This double element is shown 
schematically in Figure (3) (a tvi'o dimensional element is shovm only for 
convenionce) . It consists of two identically shaped element occupying the same 
physical space and having the same nodal numbers but different element numbers 
and dlfverent material identification cards (MiVTl for 3-d elements and MAT2 
for 2-d elements). Our goal is to let one of the overlapping elements form the 
Term-1 elastic contribution of Eq. (12b) and tim other form the Terra-2 visco- 
elastic dissipation contribution. In agreement with the notation of Figure (3), 
we refer to the first overlapping element as the "elastic element" and the 
second ns the "massless dissipation element". The only port that remains is to 
define the input constants on the M>\T1 (or MAT2 ) cards so that Term-1 and 
Term-2 In equations (12b) are properly formed. More specifically, the rationale 
for the selection of the MATl (or MAT2 ) constants follows from seeking out a 

set of parameters for the PplG] matrix in equation (10), (namely the Pp, and 

I.J lit 

compouGiiLs oC the e.Lastic [G]'-* matrix which are coatroliecl bv the user through 

e 

selection of the input variables on the appropriate M/\Ti cards) that will 

result in the desired [G]^ matrix In Term-2 of equation (12b)* The treatment 

is slightly different for three or two dimensional elements, consequently v^e 
treat them one at a time. 


Three DlmensionaX Viscoelastic Elements 


• Tlie elastic element contribution Is obtained by setting the following 
parameters on a Mi\TI card: 


1 ) 

2 ) 


set mass density (UHO) - actual mass density of viscoelastic 

material 


set loss factor (GE) = 0.0 


(13) 


3) set E = p’-'(3.\^+2n''')/(X^'+p^'') 


G 


r 

y 


0 The massless dissipation element contribution, is obtained by setting the 
following parameters on another MATl card specially earmarked for this 
second overlapping element: 

1) set mass density (RUO) 0*0 (Eero to avoid double counting) 

2) set loss factor ^ (GE) X^/e (14) 
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(1-^) cont'd) 


3) set E « eR(3+2R)/(l+R) 

G = eR 
with R = 

where e Is an arbitrary parameter that is selected smalt as desired (but not 
zero). The small e parameter removes the elastic contribution of the massless 
dissipation element. A suggested value for e is that it is on the order of 
10,000 times smaller than the larger of the real components or 

The interested reader is invited to back substitute the above values for E, 
G, and GE (e.g. rip) into the matrix formed by NASTRAN in conjunction 

hi hi G 

with equation (1) where it can be easily verified that the results reduce 
(independent of the e choice) exactly to the desired two independent parameter 

matrix [G]^ employed in Eq. (12b). The smallness of e only effects the unwanted 

elastic stiffness contribution of the massless dissipation element already 
accounted for in the "elastic element". 


Two Dimensional Elements 


• The elastic element contribution is obtained by setting the following 
parameters on a MAT2 card! 


1) set mass density (RHO) 

2) set loss factor (GE) = 

3) set Gil = 

G12 = 


= actual mass density of viscoelastic 
material 

0.0 


(15) 


G13 = 0.0 
G22 = 

G23 = C.O 
G33 = \i~ 


0 The massless dissipation element contribution is obtained by setting the 
following parameters on another MAT2 card specially earmarked for this second 
overlapping element : 


1) set mass density (RHO) ==0,0 

2) set loss factor t\^ = (GE) = X^/e 
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3) set Gll = e(l+2|iVA^) 


(16) cont'd 


G12 = e 
G13 = 0.0 
G22 = e 
G23 =0.0 
G33 = £M^/X^ 

where as with the three dimensional element, the e is an arbitrary para-- 
meter that is selected small as desired (but not zero) , See 3-d element write- 
up for a suggested e value. 

DEMONSTRATION PROBLEM 

Next we consider a demonstration problem to illustrate the implementation 
and accuracy of the two independent parameter loss factor viscoelastic elements 
described in the previous section. The problem is to determine the scattered 
and transmitted pressures for a water-submerged steel plate, (covered with a 
constant thickness layer of viscoelastic material), subject to an incident 
plane wave normal to the plate as illustrated In figure (1). The problem is 
tractable from a closed form solution point of view, consequently, an independ- 
ent check on the NASTRAN solution is available. Furthermore, experimental 
results are also available to further back up the accuracy of the physical 
representation of the viscoelastic material. 

Exact Solution 


The exact solution to this problem can initially be treated as an ordinary 
one dimensional, small deformation wave propagation problem. The effect of 
viscoelasticity can be introduced by replacing the wave speed c^ with = 

(1+iri^)^. The solution to the problem illustrated in figure (1) is carried 

out exactly like the problem given in ref. (5), page 136, except that two finite 
thick plates (rather than one) is submerged in the fluid. The origin (at x =0) 
is located at the right face of the viscoelastic layer (+x to the left). The 
back side fluid is denoted as media (4), steel plate as (3), the viscoelastic 
layer as (2) and the front side fluid as (1). The thickness of the steel plate 
is and the viscoelastic layer is ^^e following relations define the 

various waves present in the problem: 


(p.)^ . 

(P ), . B, 

si 1 


Incident wave in (1) 
Scattered wave in (1) 


(17) 
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Transmitted wave in (2) 


<»t>2 

” ^2 

i(mt-k.,x) 
e / 


“ «2 

(o»t+k2x) 

(Pt>3 

^ ^3 

^i(ut-k2(x-A2)) 


“ «3 

(wt+k^ (X-II2)) 



^i(u)t-k^(x-il2”^ 


Scattered wave in (2) 

Transmitted wave in (3) 
Scattered wave in (3) 

)) Transmitted wave in (4) 


(X7 cont'd) 


where in the fluid p is pressure and in the solid, p is the negative of the 
stress in the x direction. The k quantities are defined as: 

ki - 

"3 ■ “^=<i3 
^ ■ “/=d4 


where c^ 2 > ^^ 3 * ^^14 real dilatational wave speeds of nhe four 

materials and is the dilatational loss factor (equation ( 5 )) of the visco*- 

elastic layer. The six unknown A 29 ^2* ^ 3 * ^ 3 * ^4 determined from 

equating pressure and particle velocities at the three interfaces, thus pro- 
viding six equations to balance the six unknowns. All response variables are 
referred to the incident wave amplitude A^^, thus A^ is not considered an 
unknown in the problem. 


The quantities of interest are the scattered pressure in the incident side 
fluid (media 1 ) and the transmitted pressure in the back side fluid (media 4). 
A^ter algebraically solving for the constants we obtain 

Scatter pressure amplitude H = (^2^12 -C 1 A 22 )/ 

^^ 21 ^ 12 “‘^ 11 ^ 22 ^ ( 19 ) 

Transmitted pressure amplitude = ^ 4 /^^^ = ^^21^32^l"^32^11^2^^ ^^21^12“\l^22^ 


where 


C 


2 


= -cos(k2a2)+x r3^2Sin(k2S,2>; ^i2=P2^d2^Pl‘^dr ^34=P4^d4^Ps'^d3 
° ^ 23 ! Sxn(k2S-2)-’^23’^12 ’ ’^23^'^3‘^d3'^^2'^d2’ '^d 2 ^‘^d 2 ' 
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A ^2 * ([l+r 3 ^][cos( 2 k 3 )l 3 )+i Sin( 2 k 3 fc 3 ) ] / [ 1 ^ 3 ^] )-l 
A 21 = -r23l Sin(k2£2)-^-23^12 
A^j^ = cos(k2^.2)+i Sin(lc2^2^ 

A 22 “ 1 + tl+r3^Hcos(21c3Jl3)+l Sin(2k3A3)]/[l-r3^] 

A 32 = -2r3^[co3(k3ft3) + i Sin (k 3 Jl 3 )]/ [ 1 -^ 3 ^] 

Finite Element Solution 

The finite element representation of the figure (1) problem is shown in 
figure (2). The procedure for representing the infinite domain of fluid on the 
front and back side of the submerged plate is given in detail in ref. ( 6 ). 
Briefly, it consists of using a plane wave boundary condition at the mesh 
termination of the form p=pc^u where u is the normal particle velocity and p is 

the total pressure at the mesh termination. The nodes along the outer bound- 
aries are constrained to move only in the direction of wave propagation. The 
figure ( 2 ) sketch is drawn to scale and represents the actual number of elements 
used in the problem. All elements employed in the model are comprised of GQDMEM 
quadralateral elements. The viscoelastic zone is made up from the overlapping 
double elements described earlier in the paper (e.g. a typical viscoelastic 
element is Illustrated in figure (3)), 


Comparative Results 

The demonstration problem (both analytical and finite element) was evaluated 
with the following set of physical constants 

Table 1 - DEMONSTRATION PROBLEM PHYSICAL CONSTANTS 


MATERIAL 

psi 

r 

P 

psi 

psi 

1 

w 

psi 

P 

lb-sec^ /in^ 

Water 

345,600. 

0.0 

0.0 

0.0 

.000096 

Steel 

17,307,000. 

11,538,000. 

0.0 

0.0 

.000735 

Viscoelastic 

Material 

86,703. 

115,9 

41,736.8 

11.6 

.0003599 


The viscoelastic constants were evaluated by W, Madaigosky at NSWC. For the 
NASTRAN computer run, the water and steel plate properties were entered on a 
MAT2 card. The elastic entries Gll, Gl2, etc. correspond to the column and 
row data in the first, second and sixth (columns and rows) of the [G]^ matrix 
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of equation (11). The viscoelastic data were entered on two separate MAT2 
cards corresponding to the procedure described by equations (15) and (16). An 
E parameter oC e = A. 16 was used in the actual run. The non-dimensional results 
are shown in Table 2. 

The significance of the added viscoelastic layer on the scattered and 
transmitted pressure is seen by comparing the results of Tables 2 and 3, where- 
in the viscoelastic layer has the effect of absorbing a significant amount of 
the incident energy. The accuracy of the NASTRAN results for both the case 
with the viscoelastic layer and without the layer are quite good. As expected, 
the finite element solution accuracy falls off as the incident frequency 
increases. As pointed out by ref. (7), for problems of this type, at least 8 
elements per wave length are needed to accurately model elastic waves in the 
media. Note in Table 2 Chat as the incident frequency approaches the 8 
elements per wave length limit, the accuracy is starting to deteriorate. The 
8 element per wave length limit suggested by ref. (7) was in the absence of 
structural damping; perhaps the results presented here suggests that more than 
8 elements per x^ave length are required for structural damping (e.g. 12 
elements per wave length) . 


CONCLUDING REMARK? 

The procedure outlined here provides the NASTRAN user with an expanded 
structural damping capability, thus permitting the user to specify two independ- 
ent loss factors rjg, via the construction of the "overlapping double visco- 
elastic element" process described in this papei . The demonstration problem 
fc'hoxsTs good accuracy of the procedure relative ti. the exact solution for the 
same problem. It is acknowledged that there is some added solution time due to 
the added calculation time for the formation of the stiffness of the second 
overlapping element; however, this is a relatively insignificant amount in 
comparison to the overall solution time of the problem. The double element 
uses the same node numbers for both the "elastic element" and the superimposed 
"massless dissipation element", consequently, the matrix size or bandwidth 
properties are not affected at all. 
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Table 2 - COMPAEATIVE SOLUTION RESULTS 

(with viscoelastic layer present) 


N and imens ional 
frequency ; 

aiJ-s/Cdi 

NASTRAN 

Scattered 

Pressure 

ratio 

Exact 
Scattered 
Pressure ratio 
(equation (19)) 

NASTRAN 

transmitted 

pressure 

ratio 

Exact 

transmitted 

pressure 

ratio 

Exp er imen t al 

Scattered 

pressure 

Elements per 
wave length 
in elastomer 

.6545 

.401 

,398 

.273 

.273 

_ 

158, 

3.272 

.191 

,190 

.028 

.028 

- 

31,6 

6.545 

.108 

.104 

.006 

.007 

— 

15.8 

9.817 

,122 

,113 

.003 

.003 

.12 

10-5 

13.089 

.133 

.113 

.002 

.003 

.12 

7.9 


Note that ^ 3/^2 ^ 2,2727 (ratio of steel plate-to-viscoelastic layer thickness) 


Table 3 - COMPARATIVE SOLUTION RESULTS 

(without viscoelastic layer present) 


Nondimensional 

frequency; 

cuis/Cdi 

NASTRAN 

Scattered 

Pressure 

ratio 

Exact 
Scattered 
Pressure ratio 
(equation (19)) 

NASTRAN 

transmitted 

pressure 

ratio 

Exact 

transmitted 

pressure 

ratio 

Experimental 

Scattered 

pressure 

- 654 

.930 

.928 

.374 

.373 


3.27 

,995 

.996 

.090 

.090 

- 

6.54 

.998 

.998 

.067 

-067 

- 

9.82 

.996 

.994 

.113 

.113 

- 

13.09 

.957 

.952 

.306 

.307 

- 
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FiRure 1. Submerged Plate With Viscoelastic Layer 



Figure 2. One Dlaiensional Finite Element Model 



Figure 3. Overlapping Double Viscoelastic Element 
(Quadralateral Elanent Example) 
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